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Abstract 

We present a program for the numerical evaluation of multi-dimensional poly- 
nomial parameter integrals. Singularities regulated by dimensional regularisation 
are extracted using iterated sector decomposition. The program evaluates the co- 
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and Computing. 



Nature of problem: 

Extraction of ultraviolet and infrared singularities from parametric integrals ap- 
pearing in higher order perturbative calculations in gauge theories, e.g. multi-loop 
Feynman integrals, Wilson loops, phase space integrals. 

Solution method: 

Algebraic extraction of singularities in dimensional regularisation using iterated 
sector decomposition. This leads to a Laurent series in the dimensional regularisa- 
tion parameter e, where the coefficients are finite integrals over the unit-hypercube. 
Those integrals are evaluated numerically by Monte Carlo integration. 
Restrictions: Depending on the complexity of the problem, limited by memory and 
CPU time. Multi-scale integrals can only be evaluated at Euclidean points. 

Running time: 

Between a few minutes and several days, depending on the complexity of the prob- 
lem. 
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LONG WRITE-UP 



1 Introduction 

Sector decomposition is an algorithmic method to isolate divergences from 
parameter integrals as they occur for instance in perturbative quantum field 
theory. Originally it was devised by Hepp [I] in the context of the the proof of 
the BPHZ theorem in order to disentangle overlapping ultraviolet singularities. 
Similar ideas, applied to the subtraction of infrared divergences, can be found 
e.g. in [2]. It was employed later to extract logarithmic mass singularities from 
massive multi-scale integrals in the high energy limit at two loops [3"f4"] . 

In [5] , the concept of sector decomposition was elaborated to a general algo- 
rithm in the context of dimensional regularisation, allowing the isolation of ul- 
traviolet as well as infrared singularities from Feynman parameter integrals in 
an automated way. First applications of this algorithm were the numerical eval- 
uation of two-loop box diagrams at certain Euclidean points, see e.g. [5H6fT] . 
More recently, the method has been used to numerically check a number of an- 

most of them produced by either the public program FIESTA [22"f2"3"] or the 
code which is described in the present article. Further references about recent 
applications of sector decomposition to multi-loop calculations can be found 
in [25123] . 

Sector decomposition also has been combined with other methods for a numeri- 
cal calculation of loop amplitudes, first on a diagrammatic level in Refs. [232S], 
later for whole amplitudes in Refs. [27|2 8.29.30j. The latter approaches contain 
a combination of sector decomposition and contour deformation [31f32f33f34f35 
which allows one to integrate the Feynman parameter representation of an am- 
plitude numerically in the physical region. 

As phase space integrals in D dimensions can be written as dimensionally reg- 
ularised parameter integrals, sector decomposition can also serve to factorise 
entangled singularity structures in the case of soft and collinear real radiation. 
This idea was first presented in [36] and was subsequently applied to calculate 
all master four-particle phase space integrals where up to two particles in the fi- 
nal state can become soft and/or collinear [37]. Shortly after, this approach has 
been extended to be applicable to exclusive final states as well by expressing 
the functions produced by sector decomposition in terms of distributions [3B] • 
Further elaboration on this approach [39PU] has lead to differential NNLO re- 
sults for a number of processes [3TP2|4 3,44,45 ,46,47,48 ,49]. The combination 
of the Frixione-Kunszt-Signer subtraction scheme [50] for soft and collinear 
real radiation with the decomposition into sectors to treat real radiation at 
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NNLO, as proposed recently in [51], is also promising with regards to the 
reduction of the number of functions produced by the decomposition. A com- 
bination of sector decomposition with non-linear variable transformations as 
proposed in [52] can also serve to reduce considerably the number of functions 
to integrate, but is less straightforward to automate completely. 

To date, the method of sector decomposition has been applied successfully to 
a considerable number of higher order calculations, for a review we refer to 
[52121] • Here we will concentrate on the method of sector decomposition from 
a programming point of view. 

Despite its success in practical applications, for quite some time there was 
no formal proof for the existence of a strategy for the iterated sector decom- 
position such that the iteration is always guaranteed to terminate. This gap 
has been filled in Ref. [S3], by mapping the problem to Hironaka's Polyhedra 
game [55] and offering three strategies which are proven to terminate. Bogner 
and Weinzierl also implemented the algorithm in a public computer program 
for iterated sector decomposition written in C++ [S3]. A Mathematica inter- 
face to this program, which also allows the calculation of contracted tensor 
integrals, has recently been published in [56J. 

A different strategy guaranteed to terminate, leading to less subsectors than 
the strategies of Ref. [S3], was given by A.V. Smirnov and M. Tentyukov, who 
implemented the algorithm in the public program FIESTA [22]. Based on a 
detailed analysis of Hepp and Speer sectors in Ref. [57], an alternative strategy, 
which is based on Speer sectors, has been implemented in FIESTA2 [23]- As 
the latter strategy also uses information on the topology of the graph, it can 
perform the decomposition more efficiently in certain cases. 

Another group has implemented [SB] the sector decomposition algorithm in 
FORM jSH] . Mapping sector decomposition to convex geometry and using algo- 
rithms in computational geometry lead to a guaranteed terminating strategy 
which seems to be optimal with regards to the number of produced subsec- 
tors [60l6Tj . 

Compared to the existing packages, the present program for sector decomposi- 
tion has several new features, the main ones being a sophisticated treatment of 
(potential) numerical instabilities and the possibility to apply the program to 
more general functions, like e.g. parameter integrals occurring in phase space 
integrals over real radiation matrix elements, or functions including symbolic 
parameters. Half-integer powers of polynomial functions are also allowed. 

The structure of the article is as follows. In section [21 we briefly review the 
theoretical framework. Then we give an overview of the program structure 
and explain the individual components. The installation and usage of the pro- 
gram are described in section [31 followed by a number of examples illustrating 
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different aspects of the program in section [51 The appendix contains more 
details about timings and phase space parametrisations as well as a section 
with further information for the user. 



2 Theoretical background 



2.1 Feynman integrals 



A general Feynman graph in D dimensions at L loops with N propa- 

gators and R loop momenta in the numerator, where the propagators can have 
arbitrary, not necessarily integer powers Uj, has the following representation 
in momentum space: 



IlP?{ik},{p},7n}) 
.7=1 

4-D 



d% = d D ki , Pj({k}, {p}, m)) = q)- m) + i5 , (1) 

Z7T 2 



where the qj are linear combinations of external momenta pi and loop momenta 
h\. Introducing Feynman parameters leads to 



V(N \ °r N N r 

GZ:C = ^r^r, J n ^ *T *U - £ *«) / dV • • • 



K il • • • 



L L 

Y, kj M ij k j -2Y / kj -Qj + J + iS 



(2) 



where iV^ = X)jLi ^, M is a L x I matrix containing Feynman parameters, 
Q is an L-dimensional vector composed of external momenta and Feynman 
parameters, and J contains kinematic invariants and Feynman parameters. 

To perform the integration over the loop momenta ki, we perform the following 
shift in order to obtain a quadratic form for the term in square brackets in 
eq. ©: 



(3) 
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After momentum integration one obtains 
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where 



J r (f) = det(M) 
W(f) = det(M) 



Y,QjM- l 1 Q l -J 

.3,1=1 

l=Uv 



iS 



(5) 



and L-R/2J denotes the nearest integer less or equal to R/2. The expression 
[(M -1 (g)g)' m ) I (• R -2m)jri,...,r fl s ^ anc [ s f or sum over a \\ different combinations 
of R double- indices distributed to m metric tensors and (R — 2m) vectors I. 
The double indices Tj = (I, fii(l)) , I G {!,..., L}, i G {1, . . . , R} denote the 



;th 



Lorentz index, belonging to the I loop momentum. 



As can be seen from Eq. (jl]), the difference between scalar (R = 0) and 
tensor (R > 0) integrals, once the Lorentz structure is extracted, is given 
by the fact that there are additional polynomials of Feynman parameters 
in the numerator. These polynomials can simply be included into the sector 
decomposition procedure, thus treating contracted tensor integrals directly 
without reduction to scalar integrals. 

The functions U and J 7 also can be constructed from the topology of the 
corresponding Feynman graph. For more details we refer to j62 63 5 3f2l] . 

U is a positive semi-definite function. Its vanishing is related to the UV sub- 
divergences of the graph. Overall UV divergences, if present, will always be 
contained in the prefactor T(N U — m — LD/2). In the region where all in- 
variants formed from external momenta are negative, which we will call the 
Euclidean region in the following, J 7 is also a positive semi-definite function 
of the Feynman parameters Xj. Its vanishing does not necessarily lead to an 
IR singularity. Only if some of the invariants are zero, for example if some 
of the external momenta are light-like, the vanishing of J 7 may induce an IR 
divergence. Thus it depends on the kinematics and not only on the topology 
(like in the UV case) whether a zero of J 7 leads to a divergence or not. The 
necessary (but not sufficient) conditions for an IR divergence are given by the 
Landau equations [64 65.66], which, in parameter space, simply mean that the 
necessary condition J 7 = for an IR divergence can only be fulfilled if some 
of the parameters Xi go to zero, provided that all kinematic invariants formed 
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by external momenta are negative. 



For a diagram with massless propagators, none of the Feynman parameters 
occurs quadratically in the function J 7 = J-o . If massive internal lines are 

N 

present, J 7 gets an additional term J-{x) = J r (x)+U(x) Xjm^. If the power 

3=1 

of the Feynman parameters in the polynomial forming J 7 is larger than one for 
at least two different parameters, initially or at a later stage in the iterated 
decomposition, an infinite recursion can occur. This happens in the example 
given in section 15.1.21 if the "naive" decomposition strategy is employed. We 
have implemented a heuristic procedure to change to a different decomposition 
strategy only in cases where at least two Feynman parameters occur quadrati- 
cally, which lead to a terminating algorithm without producing a large number 
of subsectors. We do not claim that this procedure is guaranteed to terminate, 
but it proved useful for practical purposes. 



2.2 General parameter integrals 



The program can also deal with parameter integrals which are more general 
than the ones related to multi-loop integrals. The only restrictions are that 
the integration domain should be the unit hypercube, and the singularities 
should be only endpoint singularities, i.e. should be located at zero or one. 
The general form of the integrals is 



where Pi(x, {a}) are polynomial functions of the parameters Xj, which can 
also contain some symbolic constants {«}. The user can leave the parameters 
{a} symbolic during the decomposition, specifying numerical values only for 
the numerical integration step. This way the decomposition and subtraction 
steps do not have to be redone if the values for the constants are changed. 
The Vi are powers of the form z/j = a,; + b^e (with a, such that the integral is 
convergent; note that half integer powers are also possible). We would like to 
point out that most phase space integrals in D dimensions over real radiation 
matrix elements can also be remapped to functions of the type (JS)). Examples 
are given in section 




(6) 
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2.3 Iterated sector decomposition 



Here we will review the basic algorithm of iterated sector decomposition only 
briefly. For more details we refer to 021] • 

Our starting point is a function of the form of Eq. (jH]). Loop integrals in the 
form of Eq. (j3J), with open Lorentz indices contracted by external momenta 
or metric tensors, are a special case thereof, distinguished by the presence 
of a 5-distribution which constrains the sum of the integration parameters. 
Therefore loop integrals are treated somewhat differently than the more gen- 
eral functions, meaning that the 5-constraint is integrated out in a special way 
leading to the so-called primary sectors. 

I. Generation of primary sectors (loop integrals only) 

We split the integration domain into N parts and eliminate the 5-distribution 
in such a way that the remaining integrations are from to 1. To this aim we 
decompose the integration range into N sectors, where in each sector /, x\ is 
largest: 



The integral is now split into N domains corresponding to N integrals Gi from 
which we extract a common factor: G = ( — l) N,/ r(N u — LD/2)J2iLiGi. The 
N generated sectors will be called primary sectors in the following. In the 
integrals Gi we substitute 



and then integrate out x\ using the 5-constraint. As U, T are homogeneous of 
degree L,L + 1, respectively, and xi factorises completely, we have U(x) — > 
Uiit) Xi and T{x) — > ) and thus, using/ dxi/xi 5(1— Xi(l+J2k=i x k)) = 
1, we obtain 




(7) 



XiXj for j < I 
Xi for j = I 
XtXj-i for j > I 
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/ = l,...,iV. 



(9) 



8 



Note that the singular behaviour leading to 1/e-poles still comes from regions 
where a set of parameters {xj} goes to zero. 



77. Extraction of the singular factors 

For functions of the form Eq. ([6]), we first have to determine which of the 
Feynman parameters generate singularities at Xj — 1 , and which ones can lead 
to singularities at zero and one. The parameters Xj for which a denominator 
vanishes at Xj = 1 but not at Xj = should be remapped by the transformation 
Xj — > 1 — Xj. If the integrand can become singular at both endpoints of the 
integration range for a parameter Xj, we split the integration range at 1/2: 
After the split 

\ 

J dxj = J dxj + J dxj (10) 



(b) 



and the substitution Xj = Zj/2 in (a) and Xj = 1 — Zj/2 in (b), all endpoint 
singularities occur at Zj — > only. This splitting is done automatically by the 
program; the user only has to define which variables should be split. 

At this stage our starting point is a parametric integral where the integrand 
vanishes if some of the integration parameters go to zero. Our aim is to fac- 
torise the singularities, i.e. extract them in terms of overall factors of type 
x^^ 6 , aj < —1. We proceed as follows. 

1. Determine a minimal set of parameters, say S = {x^, . . . , xp r }, such that 
at least one of the functions Pi(x, {a}) vanishes if the parameters of S are 
set to zero. 

2. The corresponding integration range is an r-cube which is decomposed into 
r subsectors by decomposing unity according to 

f[ 9(1 - x Pj > 0) 6{x P] ) = t]] 9(x Ph - x Pj > 0) 0(x Pi ) . (11) 

9=1 k=l 3=1 

3. Remap the variables to the unit hypercube in each new subsector by the 
substitution 

( 12 ) 

X/3 k for j = k . 
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This gives a Jacobian factor of i» t . By construction xp k factorises from at 
least one of the functions Pi(x, {a}). 

For each subsector the above steps have to be repeated as long as a set S 
can be found such that one of the rescaled functions Pi(x, {a}) vanishes if the 
elements of S are set to zero. This way new subsectors are created in each 
subsector of the previous iteration, resulting in a tree-like structure after a 
certain number of iterations. The iteration stops if the functions Pi(x, {«}) 
contain a constant term, i.e. if they are of the form 



Pi(x, {a}) = a + Qi(x ,{a}) , (13) 

where Qi(x, {a}) are polynomials in the variables xj, and a is a constant, i.e. 
lim^o Qi{%, {«}) is nonzero. 



The resulting subsector integrals have the general form 



N 

cii 4-bje 



1 = J [u a*, ] n Pi& wr • (14) 



i=l 



The singular behaviour of the integrand now can be read off directly from the 
exponents aj, bj for a given subsector integral. 



III. Subtraction of the poles 

For a particular Xj, the integrand after the factorisation described above, is of 
the form 



i 

Ij = J dxj .r": X(x j , {x^j}, e) . (15) 
o 

If a,j > —1, no subtraction is needed and one can go to the next variable Xj + \. 
If Oj < —1, one expands I(xj, {xi^j}^ into a Taylor series around Xj = 0. 
Subtracting the Taylor series (to ordejjjp for \aj\ — p+ 1) and adding it back 
in integrated form, we obtain a part where the poles are subtracted and a 
part exhibiting 1/e poles times a function depending only on the remaining 
integration parameters. 



1 To account for half-integer exponents, e.g. dj = —3/2, we use U a jl_K denoting the 
nearest integer less or equal to \aj\. 
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R(x,e)=l(x,e)- Ua j: ' xf (0, {x^}, e) § . (16) 



For a, = —1, expanding the above expression in e is equivalent to an expansion 
in "plus distributions" [6"Tf38] 



n=0 



nl 



\n n (x) 



X 



where 
i 



dxf(x) [g(x)/x] 



X 



(17) 



with the integrations over the terms containing 5(x) already carried out. 



After having done the subtractions for each Xj, all poles are extracted, such 
that the resulting expression can be expanded in e. This defines a Laurent 



series m e 



/= E C n e n + 0(e r+1 ) , (18) 

n=-LP 

where the coefficients are finite parameter integrals of dimension (N — 1 — \ n\) 
for n < and of dimension (N — 1) for n > 0. LP denotes the leading pole, 
which can be at most 2L for an L-loop integral. The finite coefficient functions 
can be integrated by Monte Carlo integration if the Mandelstam invariants in 
T respectively the numerical constants in a general integrand have been chosen 
such that the integrand does not vanish in the integration domain. 



Improving the numerical stability 

For dj — — 1 in eq. f|T6|) . the singularity is of logarithmic nature, i.e. ~ log(A) 
if a lower cutoff A for the parameter integral was used. In renormalisable 
gauge theories, linear (a, = —2) or even higher (a, < —2) poles should not 
occur. However, they can occur at intermediate stages of a calculation, and as 
they are formally regulated by dimensional regularisation, a method has been 
worked out for the program to be able to deal with higher than logarithmic 
singularities efficiently. This method relies on integration by parts (IBP) in a 
way which aims at maximal numerical stability. 
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Let us consider an integrand after subtraction, as the one in eq. (flEl) . and 
focus on only one variable, Xj = x, and define 



I(R, a,b) = J dx x a+be R(x, e) 



(19) 



where, in the case a < 0, R(x, e) is 0(x a ) as x — > by construction. Integra- 
tion by parts gives 



I BP (R,a,b) 



a + 1 + be 
1 

a + 1 + be 



n+l+be 



R(x,e) 



dxx a+l+be R'(x 



{R(l,e)-I(R',a + l,b)} . 



(20) 
(21) 



As x a+1+be R(x, e) vanishes as x — > 0, the term ~ -R(0, e) in the square bracket 
in eq. (12"01 is zero. 

Also notice that R'(x) ~ x~ a ~ x as x — > 0, so I(R', a + 1, b) can be treated in 
the same way. 



For the case a = — 1 we thus obtain, expanding in e 



1 r 00 ffcfl" 

/Bp(i2,-1,6) = ~ dxR'(x) J2^ln n (x). (22) 

66 q n=l n - 

For a < —1, we can use eq. (I2T1) to iterate this procedure until we reach 
a >= -1. 

This method is certainly beneficial in the case of numerical instabilities coming 
from terms of the form [f(x) — f(0)]/x or [f(x) — /(0) — x f'(0)]/x 2 in R(x, e), 
leading to differences of large numbers for x — > 0. The IBP method raises the 
powers of x, trading them for additional logarithms in the integrand, which 
can be integrated more easily by the Monte Carlo program. Note also that the 
whole procedure is linear in R(x, e), so it allows one to split R(x, e) into smaller 
functions which can be dealt with more easily by the numerical integrator. 



Error treatment 

We usually integrate the sum of a small set of functions stemming from the 
e-expansion of a certain pole structure individually, and afterwards sum all the 
individual results contributing to a certain pole coefficient. The set of functions 
to sum before integration is defined by the size of the individual functions: the 
functions are summed until their sum reaches about two Megabytes. 
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The errors are calculated by adding the Monte Carlo errors of the individual 
integrations in quadrature. It is possible that there are large cancellations 
between different functions contributing to the same pole coefficient. In such 
cases it is better to first sum all coefficient functions and then integrate. The 
program offers the possibility to do so as an option which can be specified in 
the input parameter file. 

The user should be aware that for complicated functions containing many 
subtractions, the Monte Carlo error estimate is not quite appropriate: it is 
calculated on a purely statistical basis, scaling like l/y/N if N is the num- 
ber of sampling points. However, this is only a reliable error estimate under 
the assumption that the sampling has mapped all the important features of 
the function (i.e. all peaks) sufficiently precisely, and strictly is only valid for 
square integrable functions. If the function is not square integrable (but in- 
tegrable), the Monte Carlo estimate for the integral will still converge to the 
true value, but the error estimate will become unreliable. For more involved 
integrals, we are faced with functions which have gone through numerous de- 
compositions and subtractions, such that their shape in the unit hypercube is 
quite complicated, and therefore the naive Monte Carlo error estimate tends 
to underestimate the "true" error. 

Often the main source of underestimated errors in the final result is the fact 
that there are a large number of integrations to sum, and so adding the errors 
in quadrature would only give a truly appropriate error estimate if there were 
no systematic errors in the numerical integration. 

We should also note that converting the result to extract a different e-dependent 
prefactor may lead to cancellations between different contributions to a certain 
pole coefficient such that the error estimate may be too optimistic in these 
cases. 



3 Structure of the program 

The program consists of two parts, an algebraic part and a numerical part. 
The algebraic part uses code written in Mathematica [68] and does the de- 
composition into sectors, the subtraction of the singularities, the expansion in 
e and the generation of the files necessary for the numerical integration. In the 
numerical part, Fortran functions forming the coefficient of each term in the 
Laurent series in e are integrated using the Monte Carlo integration program 
BASES, version 5.1 [69J, or one of the routines from the CUBA library, ver- 
sion 2.1 [70]. The different subtasks are handled by perl scripts. The directory 
structure of the program is shown in Fig. Q], while the flowchart in Fig. [2] shows 
the basic flow of input /output streams. 
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SecDec 



basesvb.X 



Cuba 2.1 



loop 



general 



doc 



perlsrc 



demos 



perlsrc 



demos 



util 



deco 



subexp 



util 



deco 



subexp 



Fig. 1. Directory structure of the SecDec program. 

The directories loop and general have the same global structure, only some 
of the individual files are specific to loops or to more general parametric func- 
tions. The directories contain a number of perl scripts steering the decompo- 
sition and the numerical integration. The scripts use perl modules contained 
in the subdirectory perlsrc. 

The Mathematica source files are located in the subdirectories src/deco: files 
used for the decomposition, src/subexp: files used for the pole subtraction and 
expansion in e, src/util: miscellaneous useful functions. For the translation 
of the Mathematica expressions to Fortran77 functions we use the package 
Format .m [71]. The subdirectories basesv5.1 and Cuba-2.1 contain the li- 
braries for the numerical integration, taken from [69] and [70], respectively. 
The documentation, created by robodoc [72] is contained in the subdirectory 
doc. It contains an index to look up documentation of the source code in html 
format by loading masterindex.html into a browser. 

The intermediate files and the results will be stored in a subdirectory of the 
working directory whose name mysubdir can be specified by the user (first 
entry in param. input, leaving this blank is a valid option). A subdirectory of 
mysubdir with the name of the graph, respectively integral to calculate will be 
created by default. If the user would like to store the files in a directory which 
is not the subdirectory of the working directory, for example in /scratch, he 
can do this by specifying the full path in the second entry in param. input. 
An example of a directory structure created by running the examples NPbox, 
QED, ggttl, A61, a user-defined 3-loop example, and a 4-loop example to be 
written to the scratch disk is given in Fig. |3j 

The directory created for each graph will contain subdiretories according to 
the pole structure of the graph. The labelling for the pole structure is of the 
form e.g. 210h0, denoting 2 logarithmic poles, no linear and no higher poles. 
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directory loop or general 



subdirectory, e.g. graph 



param.input 
Template, m 



./launch 



makeFU.pl 
decompose.pl 



subexp.pl 





Mathematica output 
graph * .out 


» 


subtractions 




Fortran functions : 
polestructurej * ./ 

executables 


> 

expansion 

compilation 


> 


launch integration 


> 




epstothe [i] / point [i] .out 




> 


' results.pl 




graph [point] full. res 



BASES 
CUBA 



Fig. 2. Flowchart showing the main steps the program performs to produce the result 
files. In each of the subdirectories loop or general, the file Template .m can be used 
to define the integrand. The produced files are written to a subdirectory created 
according to the settings given in param.input. By default, a subdirectory with 
the name of the graph or integrand is created to store the produced functions. This 
directory will contain subdirectories according to the pole structure of the integrand. 
The perl scripts (extension .pi) are steering the various steps to be performed by 
the program. 

We should point out that this labelling does not necessarily correspond to the 
final pole structure of the integral. It is merely for book-keeping purposes, and 
is based on the counting of the powers of the factorised integration variables. 
In more detail, if %\ variables have power —1, i 2 variables have a power —2 < 
i 2 < —1 and is variables have a power < —2, the labelling will be iilizhis, 
even though the non-logarithmic poles will disappear upon e-expansion. In 
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/scratch 



2loop 



3loop 



-Hoop 



NPbox 




QED 




ggttl 




A61 







Fig. 3. Example for a directory structure created by running the loop demo programs 
NPbox, QED, ggttl, A61. A four-loop example denned by the user to be written to 
the scratch disk is also shown. 



QED 





2/0ft0 
















epstothe — 2 


epstothe — 1 


epstothe 



UOhO 



epstothe — 1 



epstothe 



0/0/i0 



epstothe 



Fig. 4. Example for a directory tree corresponding to the pole structure of the graph 
QED contained in the demo programs. 

particular, for half-integer powers, the labelling does not correspond to "true" 
poles, but rather to terms which can be cast into functions like T(— 3/2 — e), 
which are well-defined in the context of dimensional regularisation, where e 
can be regarded as an arbitrary (complex) parameter. Note also that in the 
case of a prefactor containing 1/e poles multiplying the parameter integral, the 
poles which are nagged up at this stage of the program will only correspond 
to the poles read off from the integration parameters. In any case, the final 
result will be given to the order specified by the user in param. input. 
Each of these "polestructure" directories contains further subdirectories where 
the files for a particular power in epsilon are stored. An example is given in 
Fig.H 



The user only has to edit the following two files: 



• param. input: (text file) 

specification of paths, type of integrand, order in e, output format, param- 
eters for numerical integration, further options 

• Template. m: (Mathematica syntax) 

• for loop integrals: specification of loop momenta, propagators; optionally 
numerator, non-standard propagator powers 
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• for general functions: specification of integration variables, integrand, vari- 
ables to be split 

To give a specific example rather than empty templates, the files param . input 
and Template. m in the loop subdirectory contain the setup for example 1, 
the non-planar massless on-shell two-loop box diagram, while those in the 
general directory contain the setup for example 6, a hypergeometric function 
of type 5 F 4 . 

Apart from these default parameter /temp late files, the program comes with 
example input and template files in the subdiretories loop/demos respectively 
general/demos, described in detail in section |5j 

The user can choose the numerical integration routine and the settings for the 
different integrators contained in the Cuba library in the file param. input. 
The compilation of the chosen integration routine with the corresponding set- 
tings will be done automatically by the program. 



4 Installation and usage 



Installation 



The program can be downloaded from http: //projects .hepforge . org/secdec/ 



Installation is done by unpacking the tar archive, using the command tar xzvf 
SecDec.tar.gz. This will create a directory called SecDec with the subdirecto- 
ries as described above. Change to the SecDec directory and run ./install. 

Prerequisites are Mathematica, version 6 or above, perl (installed by default 
on most Unix/Linux systems) and a Fortran compiler (e.g. gfortran, ifort). 
The install script only checks if Mathematica and perl are installed on the 
system and inserts the corresponding path into the perl scripts. The install 
script does not test the existence of a Fortran compiler because the compiler 
should be specified by the user in param. input. If no compiler is specified, it 
defaults to gfortran. 



Usage 

(1) Change to the subdirectory loop or general, depending on whether you 
would like to calculate a loop integral or a more general parameter inte- 
gral. 
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(2) Copy the files param. input and Template .m to create your own param- 
eter and template files myparamf ile, mytemplatef ile. 

(3) Set the desired parameters in myparamf ile and define the propagators 
etc. in mytemplatef ile. 

(4) Execute the command ./launch -p myparamfile -t mytemplatefile in the 
shell. 

If you omit the option -p myparamfile, the file param . input will be taken 
as default. Likewise, if you omit the option -t mytemplatefile, the file 
Template. m will be taken as default. If your files myparamfile, mytem- 
platefile are in a different directory, say, myworkingdir, use the option 
-d myworkingdir, i.e. the full command then looks like ./launch -d my- 
workingdir -p myparamfile -t mytemplatefile, executed from the directory 
SecDec/loop or SecDec/general. 

Alternatively, you can call the launch script from any directory if you 
prepend the path to the launch script, i.e. the command 
pathdo -launch /launch -p myparamfile -t mytemplatefile executed from 
myworkingdir would run the program in the same way. pathdo -launch 
can be either the full or relative path for SecDec/loop or SecDec/general. 

The ./launch command will launch the following perl scripts: 

• makeFU.pl: (only for loop integrals) constructs the integrand functions 
J 7 , U and the numerator function from the propagators and indices given 
in Template .m. 

• decompose.pl: launches the iterated sector decomposition 

• subexp . pi: launches the subtractions and epsilon-expansions and writes 
the Fortran functions. Depending on the "exe-flag" specified in the pa- 
rameter file (see below for a detailed explanation of the flag), this script 
also launches the compilation and the numerical integrations. 

(5) Collect the results. Depending on whether you have used a single machine 
or submitted the jobs to a cluster, the following actions will be performed: 

• If the calculations are done sequentially on a single machine, the results 
will be collected automatically (via results.pl called by launch). The 
output file will be displayed with your specified text editor. The re- 
sults are also saved to the files [graph] _ [point] epstothe* . res and 
[graph] _ [point] full . res in the subdirectory subdir/graph (loops) 
respectively subdir/integrand (general integrands) (name specified in 
param. input, where you can also specify different names for different 
numerical points). 

• If the jobs have been submitted to a cluster: when all jobs have finished, 
execute the command ./results. pi [-d myworkingdir -p myparamfile] in 
a shell from the directory SecDec/loop or SecDec/general to create 
the file containing the final results. 

If the the user needs to change the batch system settings: manually 
edit perlsrc/makejob.pm and perlsrc/launchjob.pm. This writes 
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the desired syntax to the scripts job [polestructure] in the corre- 
sponding subdir/graph or subdir/ integrand subdirectory. 

(6) After the calculation and the collection of the results is completed, you 
can use the shell command ./launchclean [graph] to remove obsolete files. 
If called with no arguments, the script only removes object files, launch 
scripts, makefiles and executables, but leaves the Fortran files created 
by Mathematica, so that different numerical points can be calculated 
without rerunning the Mathematica code. If called with the argument 
'all' (i.e. ./launchclean [graph] alt), it removes everything except the result 
files displaying the final result and the timings. 

The 'exe' flag contained in p ar am . input offers the possibility to run the 
program only up to certain intermediate stages. The flag can take values from 
to 4. The different levels are: 

exe=0: does the iterated sector decomposition and writes files containing 
lists of subsector functions (graphsec* . out) for each pole structure to the 
output subdirectory. Also writes the Mathematica files subandexpand* .m 
for each pole structure, which serve to do the symbolic subtraction, epsilon 
expansion and creation of the Fortran files. Also writes the scripts 
batch [polestructure] which serve to launch these jobs at a later stage. 

exe=l: launches the scripts batch [polestructure] . This will produce the 
Fortran functions and write them to individual subdirectories for each pole 
structure. 

exe=2: creates all the additional files needed for the numerical integration. 
exe=3: compilation is launched to make the executables. 
exe=4: the executables are run. 

If the first steps of the calculation, e.g. the decomposition or the creation of 
the Fortran functions, are already done, the following commands are available 
to continue the calculation without having to restart from scratch: 

• finishnumerics.pl [-d myworkingdir -p myparameterfile] : 

if the 'exe' flag in p ar am . input resp. myworkingdir /myparamet erf ile is 
set smaller than four, this will complete the calculation without redoing 
previous steps. 

• justnumerics.pl [-d myworkingdir -p myparameterfile]: 

if you would like to redo just the numerical integration, for example to pro- 
duce results for a different numerical point or to try out a different number 
of sampling points, iterations etc. for the Monte Carlo integration: change 
the values for the numerical point resp. the settings for the Monte Carlo 
integration and the "name of the numerical point" in the parameter file, 
and then use the command ./ justnumerics.pl [-d myworkingdir -p mypa- 
rameterfile] to redo only the numerical integrations (if the Fortran files f*.f 
have been produced already). Using this option skips the Mathematica sub- 
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traction and epsilon expansion step which can be done once and for all, 
as the variables at this stage are still symbolic. After completion of the 
numerical integrations, use the command ./results. pi [-d myworkingdir -p 
myparameterfile] to collect and display the results as above. 

The program tries to detect the path to Mathematica automatically. In case 
you get the message "path for Mathematica not automatically found" , please 
insert the path to Mathematica on your system manually for the variable 
$mathpath in the file perlsrc/mathlaunch.pl. 

We also should mention that the code starts working first on the most compli- 
cated pole structure, which takes longest. This is because in case the jobs are 
sent to a cluster, it is advantageous to first send the jobs which are expected 
to take the most time. 



5 Description of Examples 

5.1 Loop integrals 

The examples described below can be found in the subdirectory loop/demos. 



5.1.1 Example 1: Non-planar massless two-loop box 

The non-planar massless two- loop box is a non-trivial example, as the sector 
decomposition applied to the standard representation, produced by combin- 
ing all propagators simultaneously with Feynman parameters, exhibits "non- 
logarithmic poles" (i.e. exponents of Feynman parameters < —1) in the course 
of the decomposition. We should point out that, even though the program can 
deal with linear or higher poles in a completely automated way, it is often 
a good idea to investigate if the integrand can be re-parametrised such that 
poles of this type do not occur, because these poles require complicated sub- 
traction terms which slow down the calculation. Non-linear transformations 
as e.g. described in [52] can be useful in this context. Further, integrating 
out first one loop momentum, and then combine the remaining propagators 
with the obtained intermediate result using another set of Feynman parame- 
ters often leads to a representation where at least one of the parameters can 
be factorised without sector decomposition, thus speeding up the calculation 
considerably. This is demonstrated for the non-planar massless two-loop box 
in appendix 18.21 and the template to calculate the graph in this way can be 
found in SecDec/general/demos. 
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i'l 



Pi 



P2 

Fig. 5. The non-planar two-loop box, called NPbox in example 1. 

To obtain results for the non-planar massless two-loop box shown in Fig. [5] 
without doing any analytical steps, copy the file loop/param. input to a 
new parameter file, say paramNPbox. input, and specify the desired order 
in e, the numerical point and possible further options. Likewise, copy the 
file loop/Template .m to a new template file, say templateNPbox.m (this al- 
ready has been done for the examples described here, see the subdirectory 
loop/demos). Then, in the loop/demos directory, use the command ../launch 
-p paramNPbox. input -t templateNPbox.m. The file templateNPbox.m already 
has the propagators of the non-planar double box predefined. In paramNPbox . input, 
we defined the prefactor such that a factor of — r(3 + 2e) is not included in 
the numerical result: if we define 

G NP (s, t, u) = -r(3 + 2e) £ % + 0(e) , (23) 

n=0 6 

the program should yield the results for P n , given in Tabled] Note that ac- 
cording to eq. jl]), we always divide L-loop integrals by (m~) L , so this factor 
is never included in the numerical result. The decomposition produces 384 
subsectors. 



(s,t,u) 


(-1,-1,-1) 


(-1,-2,-3) 


P_ 4 


1.75006±1.3 x 10" 4 


0.41670±1.1 x 10" 4 


P-3 


-2.99969± 0.00055 


-0.9313 ±0.00067 


P-2 


-22.821 ± 0.003 


-5.8599 ± 0.0035 


P-l 


113.629 ±0.013 


42.79 ±0.02 


Po 


-395.27 ± 0.05 


-162.73±0.09 



Table 1 

Numerical results for the points (s,t,u) = (—1,-1,-1) and (—1,-2,-3) of the 
massless non-planar double box. 

The result for the graph called NPbox at the numerical point called point 
in the input file will be written to the file NPbox_ [point] full .res in the 
subdirectory 21oop/NPbox, where 21oop is a subdirectory which has been 
created by the program, using the directory name the user has specified in the 
first entry of paramNPbox. input. By default, a subdirectory with the name 
of the graph is created, but the user can also specify a completely different 
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3 6 

Fig. 6. Blue (solid) lines denote massive particles. 



directory (e.g. scratch) where the results will be written to (second entry in 
paramNPbox . input) . 

More information about the decomposition is given in the file NPboxOUT. info. 
Information about the numerical integration is contained in the files 
[point] intf ile . log in the subdirectories graph/polestructure/epstothe [i] , 
where "polestructure" is of the form e.g. 210h0, denoting 2 logarithmic poles 
and linear, higher poles. 

It should be emphasized that in par am. input, the numbers for the Mandel- 
stam invariants should be defined as the Euclidean values, so the values for 
s,t,u,p 2 should always be negative in param. input. Note also that the con- 
dition s + t + u = cannot be fulfilled numerically in the Euclidean region, 
so it should not be used in onshell={ . . . } in the template file to eliminate u 
from the function T in the case of non-planar box graphs. 



5.1.2 Example 2: Planar two-loop ladder diagram with massive on-shell legs 

The purpose of this example is to show how to deal with diagrams where 
the decomposition could run into an infinite recursion if the default strat- 
egy is applied. The rungs of the ladder are massless particles (e.g. photons), 
while the remaining lines are massive on-shell particles, depicted by solid 
(blue) lines in Fig. El To run this example, execute the command ../launch 
-p paramQED. input -t templateQED.m from the loop/demos directory. Only 
the primary sectors number one and seven are at risk of running into infinite 
recursion, therefore they are listed in the third-last item of paramQED . input 
as the ones to be decomposed by a different strategy. The results for the nu- 
merical point called point will be written to the file QED_ [point] full. res 
in the subdirectory 21oop/QED. Numerical results for some sample points are 
given in Table|2j The kinematic points are defined by the mass m and the 
Mandelstam variables s = (pi +P2) 2 ,t = (p 2 + p^) 2 ■ We extracted a prefactor 
of T(l + e) 2 . 



22 



(s,t,m) 


/no n q 1 \ 
{-\j.Z,-yj.o,l ) 


( Q /O /I /Q 1 / CC^\ 


"-2 


-l.ODlOllt I.OOXIU 


1C17 _L n nnnQ 
-z.ioi I ± U.UUUo 




r 0070 _i_n nm q 
-O.oo/o ztU.UUlo 


1 /ivm _i_n nnofi 
-1.4/U1 ±U.UUzb 


Po 


1.419 ±0.025 


30.191± 0.014 


Pi 


62.46 ± 0.18 


140.73±0.057 


P2 


284.76 ± 0.87 


450.67±0.19 



Table 2 

Numerical results up to order e 2 for the points (s,t,m) = (—0.2,-0.3,1) and 
(—3/2, —4/3, 1/5) of the two-loop ladder diagram shown in Fig.[6j An overall factor 
of r(l + e) 2 is not included in the numerical result. 
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Fig. 7. Non-planar graphs occurring in the calculation of gg — > tt at NNLO. Blue 
(solid) lines denote massive particles. 



5.1.3 Example 3: Non-planar two-loop diagrams with two massive on-shell 
legs 

This example gives results for two non-planar graphs occurring in the calcula- 
tion of gg — > tt at NNLO, shown in Fig. [71 The analytic results for these graphs 
are not yet available. Numerical results at Euclidean points can be produced by 
choosing numerical values for the invariants s,t,u,m 2 in paramggttl . input 
respectively paramggtt2 . input and then executing the command ../launch -p 
paramggttl. input -t templateggttl.m in the loop/demos directory, analogously 
for ggtt2. Results for two sample points are shown in Table |3j 



5.1.4 Example 4-' A rank one tensor two-loop box 

In order to demonstrate how to run the program for integrals with non-trivial 
numerators, we give the example of a rank one planar massless on-shell two- 
loop box, where we contract one loop momentum in the numerator by 2p%. 
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ggttl 




(-0.5,-0.4,-0.1,0.17,0.17) 


(-1.5,-0.3,-0.2,3,1) 


Po 


-38.0797±0.0027 


-0.19904±1.5 x 10~ 5 


Pi 


-263. 22± 0.015 


-0.71466±6 x 10" 5 


P2 


-936.86± 0.06 


-1.45505± 0.0002 


ggtt2 


(s,t,u,mj,ml) 


(-0.5,-0.4,-0.1,0.17,0) 


(-1.5,-0.3,-0.2,3,0) 


P-4 


-10.9159 ± 0.0006 


-0.13678±1.46 x 10~ 5 


P-3 


-43.5213 ± 0.0075 


-0.2087 ±0.00024 


P-2 


165.384 ± 0.048 


3.3417 ±0.0014 


P-l 


20.842±0.268 


-6.593±0.007 


Po 


2117.5 ± 1. 57 


20.42±0.04 



Table 3 

Numerical results for the diagrams shown in Fig. F7J The finite diagram ggttl has 
been calculated up to order e 2 . An overall factor of T(l + e) 2 is extracted. 

rd D kd D l 2p 3 -k 

(24) 

where we omitted the i5 terms in the propagators. The result for the kinematic 
sample point (s,t,u) = (—3,-2,5) is shown in TableHl Note that in this 



(s,t,u) 


(-3,-2,5) 


P-4 


-0.319449 ±1.7 x 10~ 5 


P-3 


0.46536 ±8 x 10~ 5 


P-2 


0.5848 ±0.0004 


P-l 


-3.3437 ± 0.0013 


Po 


-1.6991± 0.0035 



Table 4 

Numerical results for the point (s, t, u) = (—3, —2, 5) of the rank one two-loop ladder 
diagram given by eq. (I24p . An overall factor of T(l ± e) 2 has been extracted. 

example, we used a positive value for the Mandelstam invariant u, which 
seems to contradict the requirement to have only Euclidean values for the 
invariants. However, in this case we can do this because the function J 7 does 
not depend on u at all. The numerator does depend on u, but as a numerator 
which is not positive definite does not spoil the numerical convergence, we can 
as well choose a numerical value for u such that the relation s±t±n = 0is 
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Fig. 8. The three-loop vertex diagram A$ i with the dotted propagator raised to the 
power 1 + e. 

fulfilled. This has the advantage that it allows us to use the latter relation to 
simplify the numerator. 

5.1.5 Example 5: A three-loop vertex diagram with e-dependent propagator 
powers 

This example shows how to calculate diagrams with propagator powers dif- 
ferent from one. The results for the graph A 61 (notation of Ref. [5]), given 
in Table can be produced by running ../launch -p paramA61. input -t tem- 
plateA61.m from the loop/demos directory. 



P-3 


P-2 


P-i 


Po 


Pi 


P2 


0.16666 


1.8334 


18.123 


125.32 


889.96 


5325.3 



Table 5 

Numerical results for the diagram shown in Fig. [8]with the dotted propagator raised 
to the power 1 + e. The errors are below one percent. 

The analytical result for this diagram with general propagator powers is given 
in Ref. [8] and is also given in the file 31oop/A61/A61analytic .m to allow 
comparisons between analytical and numerical results for arbitrary propagator 
powers. 

5.2 More general polynomial functions 

The examples described below can be found in the subdirectory general/demos. 
5.2.1 Example 6: Hypergeometric functions 

As an example for "general" polynomial functions, we consider the hypergeo- 
metric functions p F p _i(ai, . . . , a p ; b±, . . . , b v -\ \ (3), using the integral represen- 
tation recursively: 
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pF p -i(ai, . . . ,a p ;b\, . . . , fe p _i; p) 
i 



'p-ij 



r(a p )r(6 p _i - Op) 

dz (1 - z )-i-«p+6p-i p _ 1 F p _ 2 (a 1 , . . . , a p _i; 61,..., 6 P _ 2 ; P) 



(25) 



T(6i 



r(a 2 )r(6i - a 2 ) 



cfe (1 — z] 



-1— a2+bi _— 1+02 



Considering 5 F 4 (ai, . . . , a 5 ; 61, . . . , 64; /3) with the values ai = e, a 2 = — e, a 3 = 
— 3e, a 4 = — 5e, a 5 = — 7e, &! = 2e, 6 2 = 4e, 63 = 6e, 64 = 8e, /3 = 0.5 we obtain 
the results shown in Table O The "analytic result" has been obtained using 
HypExp[7S|[75. 



e order 


analytic result 


numerical result 


time taken (sees) 


e° 


1 


1.0000002 ±4 x 10~ 7 


2 


e 1 


0.189532 


0.189596±0.00036 


21 


e 2 


-2.299043 


-2.306±0.011 


124 


e 3 


55.46902 


55.61 ±0.39 


248 


e 4 


-1014.39 


-1018.4±5.9 


429 



Table 6 

Results for the hypergeometric function 5F^(e, — e, — 3e, — 5e, — 7e; 2e, 4e, 6e, 8e; /3) at 
(3 = 0.5. The timings in the last column are the ones for the numerical integration. 
The time taken for decomposition, subtraction and e-expansion was 11 seconds. 

This result can be produced by typing ./launch -d demos -p param5F4 .input 
-t template5F4.nl in the subdirectory general, or by typing ../launch -p 
param5F4 .input -t template5F4-m in the subdirectory general/demos. 



The program can also deal with functions containing half integer exponents. 
Table [7] shows results for 4 F 3 with arguments a\ = — 4e, a 2 = —1/2 — e,a 3 = 
-3/2 - 2e, a 4 = 1/2 - 3e, b x = -1/2 + 2e, b 2 = -1/2 + 4e, b 3 = 1/2 + 6e. 
These results can be produced by the command ../launch -p param4F3. input 
-t template4F3.m in the subdirectory general/demos. 



e order 


analytic result 


numerical result 


time taken (seconds) 


e° 


1 


0.999997 ±1.7 x 10~ 5 


1.6 


e 1 


-4.27969 


-4.2810 ± 0.0055 


54 


e 2 


-26.6976 


-26.625 ±0.121 


90 



Table 7 

Results for the hypergeometric function 4,F^(— 4e, — 1/2 — e, — 3/2 — 2e, 1/2 — 
3e; -1/2 + 2e, -1/2 + 4e, 1/2 + 6e; /3) at f3 = 0.5. 
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5.2.2 Example 7: Phase space integrals 



Sector decomposition can be useful for the calculation of phase space integrals 
where infrared divergences are regulated dimensionally. This is particularly 
the case for double real radiation occurring in NNLO calculations involving 
massive particles, where analytic methods show their limitations. 

Here we give examples of 2 — y 3 phase space integrals, which should be con- 
sidered as part of a 2 — y n phase space written in factorised form. We choose 
particles 3 and 4 to be massless, while ps is the momentum of a massive state, 
either a single particle or a pseudo-state formed by n additional momenta pi 
in the final state, i.e. p§ = Y^?=5Pi- After all integrations have been mapped 
to the unit interval, we have integrals of the form 



J rf$ 3 =a J Hdxi - Xl )x 2 (i - x 2 )}^ ± [x 3 (i - x 3 )\ D - 

i=l 



[x 4 (1 - xj] 2 ? [1 - f3x 3 (1 - x 2 )} 2 ~ D , (26) 
P = 1 - — , C £ = —ig-j dO^-adn^ s D ^2 D ^^ . (27) 

The derivation is given in the appendix, section 18.31 
The invariants in this parametrisation are given by 

sis = ~sf3x 3 (1 - xx) 
s 23 = -s(3x 3 x 1 

s u = f3Kx 3 (l-x 2 ), K- 



S 2 A 

where 



1 - (3x 3 (1 - x 2 ) 
1 - (3(l-x 2 x 3 ) 
1 -px 3 (l -x 2 ) 
s u = ~K |t" + x A (t + - t")} = —K s u 

K |m + — x 4 (u + — w~) j = — K §24 , 



t ± =[Jx 1 (l-X 2 )± yJx 2 (l-X 1 )) (21 

n ± = ( J(l -xi)(l -x 2 ) ± ^x x x 2 



2 



We would like to point out that for the examples below, more convenient 
parametrisations, i.e. parametrisations where the variables in the denominator 
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Fig. 9. Interference of diagrams leading to factors of S35S23 in the denominator. 

factorise, and/or reflect symmetries of the squared matrix element, certainly 
do exist. However, the purpose of the examples is to illustrate that the code 
can deal with denominators which are amongst the most complicated ones 
which do occur in NNLO real radiation involving two (unresolved) massless 
particles in the final state, where they cannot always be "rotated away" by 
suitable transformations. A hybrid approach combining sector decomposition 
with convenient parametrisations/transformations is certainly the method of 
choice for real radiation at NNLO. The program can be used to evaluate the 
integrals occurring in such an approach. 



Three massless particles in the final state 

We first consider where P5 is a massless particle, i.e. the limit f3 — > 1 in 

eq. ( 1261) . If we combine the phase space with the toy matrix element 1/(335323) , 
we have singularities at X\ = 0, 22 = and £3 = 0. Such denominators come 
e.g. from the interference of diagrams as shown in Fig. |9j 

1 4 

C e Jlldxi [(l-xi)(l-z 2 )]^ ix 1 x 2 ] £ ^xt 5 a-xs) D ~ 3 
i=l 

[x 4 (1 - £4)]^ [1 - x 3 (1 - x 2 )f- D . (29) 

The term [1 — £3 (1 — X2)] 3 ~ D goes to zero for 23 — > 1, x 2 — > 0. Although £3 = 1 
does not lead to a singularity in the above example, for numerical stability 
reasons, and having in mind the presence of more complicated matrix elements 
than our toy example, it is preferable to transform this factor to an expression 
which is finite in the above limits. Splitting the £3 integration at 1/2 and 
then doing sector decomposition achieves this goal. The program will do this 
automatically if the template file contains splitlist={3}, to tell the program 
that the integration over £3 should be split at 1/2. Of course the singularities 
at Xi = will also be extracted automatically. 

Using the command . ./launch -p params23s35. input -t templates23s35.mvn the 
subdirectory general/demos, sector decomposition leads to the result given 
in Table El 



s 2 
S35S23 
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P-3 


P-2 


P-i 


^0 


-1.5705± 0.0005 


-4.3530 ± 0.0025 


1.712± 0.005 


31.040± 0.014 



Table 8 

Results for the integral given by eq. (|29p. The factor C e is not included in the 
numerical result. 

Two massless and one massive particles in the final state 



The example in this subsection illustrates the program option to exclude cer- 
tain parts of the integrand from the decomposition, even though they can 
become zero at certain values of the integration parameters. This can be use- 
ful if a particular term is known not to lead to a singularity. Note that terms 
with powers > are excluded from the decomposition by default. 

As an example we pick an integral over s 14 , where the line singularity has been 
remapped already (see appendix, section [873]) . 



/ d$ 3 lit = 2C e f d Xl dx 2 dx z dx A [xA (1 - x 4 )]^ xf~ 3 (1 - x 3 ) D ~ A (30) 

J St a J 



[1 - /3x 3 (1 - xx + x^f" X?~* x°~ 5 [(1 - Xl )(l - X 2 )(l -X X + XiX 2 )]^ 
{yj (1 - Xi)(l - X 2 ) - y/l - Xi + XiX 2 ) 2 + 4 Xa \J{1 - Xi){l - X 2 )(l - ~X\ + XiX 2 ) 



A-D 



Choosing the option "n" for "no decomposition" in the definition of the inte- 
grand for the term in square brackets [. . .} A ~ D (see templates 14. m), there will 
be no decomposition in the variables x 2 ,xa, although this term vanishes in 
the limit x 2 , xa — > 0, but this limit does not lead to a singularity. The result, 
which can be produced by ../launch -p params 14 .input -t templates 1 4 .m in 
the subdirectory general/demos, is given in Table [9l 



P-I 


^0 


-1.12635 ± 0.0003 


-8.771 ± 0.003 



Table 9 

Numerical result for the integral given by eq. (|30p for (3=0.75. The factor 2 C e is not 
included in the numerical result. 



6 Conclusions and outlook 



We have presented a program for the numerical evaluation of multi-loop inte- 
grals in Euclidean space, as well as the evaluation of more general parameter 
integrals in the context of dimensional regularisation. Singularities which lead 
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to poles in 1/e are extracted automatically using iterated sector decomposi- 
tion. The program then produces finite functions forming the coefficients of a 
Laurent series in e, which are evaluated numerically by Monte Carlo integra- 
tion. 

In the case of loop integrals, the program can deal with arbitrary tensor in- 
tegrals, where the corresponding numerator function is constructed automat- 
ically. Non-standard propagator powers, e.g. powers depending on e, are also 
supported. In the case of general polynomial functions, the program also can 
deal with cases where the integration parameters or the polynomials present 
in the integrand are raised to half-integer powers. Constants can be left sym- 
bolic at the sector decomposition stage; their values can be specified later at 
the numerical integration stage, such that the decomposition has to be done 
only once and for all while the numerical values for the constants still can be 
changed. 

The code is publicly available at |http : //pro j ects . h epf org e . org/ secdecT] 
and comes with various examples and detailed documentation. Different choices 
of numerical integration packages are possible, i.e. the Monte Carlo program 
Bases [69] and the ones from the Cuba library given in [70] . For a future ver- 
sion, we plan to include alternative tools to generate optimised functions, e.g. 
the one described in |75j . A version which includes non-linear transformations 
in an automated way, combined with the extension to integrands containing 
e.g. physical thresholds, requiring complex contour integration, is also planned. 
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8 Appendix 

8. 1 Timings for a J^-loop two-point diagram 

In order to give an idea about the timings for a complicated example which we 
ran on several processors, we give the timings for the four-loop graph shown in 
Fig.[TUJ The coefficient of each pole order is composed of a number of functions 
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Fig. 10. A four-loop two-point master integral 



which can be integrated individually, such that the time taken for the longest 
job equals the total time for a given pole order, provided that the contributing 
functions are integrated in parallel. The files to run the integration of these 
functions in parallel are created by the program automatically. The number 
of integrations to run individually depends on the size and number of the 
regular subsector functions contributing at each pole order: these functions 
are summed until their sum reaches about one Megabyte, and then integrated 
individually, preferably in parallel. 



Stage 


Time for longest job (sees) 


nb of functions 


Total time (sees) 


Sector Decomposition 


615 




2160 


Subtraction &; e-expansion 


809 




3767 


Numerical integration e" 1 


156 


5 


508 


Numerical integration e° 


422 


28 


4720 


Numerical integration e 1 


492 


28 


5946 


Numerical integration e 2 


2172 


29 


8123 



Table 10 

Timings for the diagram shown in Fig. [TUJ The time taken for the longest job equals 
the total time for a given pole order if the contributing functions are integrated in 
parallel. The number of sampling points was 500000 for each pole order. The last 
column shows the timings which would result from a calculation in series. 



The timings are listed in Table [TUJ For information we also give the timings 
which would result from a serial calculation in the third column of Table [TUJ 
The results are shown in Table [TT] 

In this calculation, the symmetry of the problem was used, so only four pri- 
mary sectors were evaluated. The corresponding multiplicities of the primary 
sectors are taken into account automatically, provided they are specified in 
par am. input. 
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Order 


Analytical result 


Numerical result 


,-i 
e 


-lU.ooyz r 


i n qvi _i_n nno 
-lU.o/ l±U.UUz 


e° 


-70.99081719 


-71.002±0.013 


e 1 


-21.663005 


-21.65±0.12 


e 2 


Unknown 


2833.79±0.92 



Table 11 

"Analytical" and numerical results for the diagram shown in Fig. [TUJ The analytical 
result can be found in |16j . 

8.2 Another representation of the non-planar two-loop box 

Here we derive a representation of the non-planar two-loop box where one 
integration parameter factorises naturally, such that it can be integrated out 
analytically, leaving a representation which can be evaluated in a completely 
automated way by the routines in SecDec/general, the evaluation being con- 
siderably faster than the one of example 15.1.11 This procedure is not limited 
to our particular example, but requires an analytical step of introducing a 
convenient parametrisation. 

The expression for the non-planar two-loop box shown in Fig. [5] is given by 



, d D kd D l 



(nrf) 2 k 2 (k + p 2 y(k- Pl y(k-iy(i + p 2 ) 2 (k-i + p 4 y(i + p 2 + p 3 y ' 

(31) 

Considering first the integration over the loop momentum I only, we have a 
one-loop box as shown in Fig. [11] with Pi — pi — k, P 2 = p 2 + k. Feynman 

Pi Ps 



Pi Pi 

Fig. 11. The "inner" box as part of the non-planar two-loop box shown in Fig.O 
parametrisation for this one-loop subgraph leads to 



J ?.7T 2 



l7l % (k - + P2 y(k - 1 + p A y{i + P2 + P3 ) 2 

r 4 4 
r(2 + 2e) / J] dxi 6{1-J2 x j) k )~ 2 ~ e 

i=l j=l 
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T(x, k) = — (pi + p 4 - /c) 2 xix 3 - (pi + P3 - ^2^4 — (A; + P2) 2 £i£2 - (pi - k) 2 x 3 x 4 
Now we substitute 



xi = t 2 (1 - t 3 ) , ^2 = *i*3 , S3 = (1 - *i) *3 (33) 
and integrate out the 5-constraint to obtain 



h = r(2 + e) | cMM£ 3 [* 3 (1 - fa)]" 



l-e 



-(P! + p 4 - /i;) 2 Ul - (Pi + P3 - /S) 2 *1*2 - (fc + P2) 2 *1<2 - (Pi - A:) 2 *1*2 



-2-e 



2r(2 + e)r(-e)r(l-e) 



1 



<iticZt 2 (34) 



r(l-2e) 

(pi + P4 - k) 2 t 2 ii - (pi + p 3 - fc) 2 tit 2 - (k + p 2 )^it 2 - (Pi - kf tih 





2 , 7 / , „ l\2j J /;„ 1 „ \2 + + /„ ?„\2 r ~ n _2 ~ e 



where we used the shorthand notation ii = 1 — tj. Now we combine the above 
expression with the remaining propagators, treating the expression in square 
brackets in eq. (134]) as a fourth propagator with power 2 + e. One can use 
the SecDec code to calculate the resulting function integrand function J-* 2 , 
although this is an easy calculation to do by hand. One obtains 



r 4 4 

3-2e 



G NP = ZL[ rn-oT I dhdt2 I II dzi z\ + %l - £ h) 

1 \ l Ze ) J J i=i 3=1 

^(z) = —S12 Z2Z3 — T Z1Z4 — P 2 Z2Z4 — P 2 2:32:4 , 
where 

T=S 23 £ 2 (1 + S 13^l(l - t 2 ) , Sy = {Vi+Vjf , 

Pi = S12 (1 " *0(1 " * 2 ) , ^ = S12 *1*2 • (35) 

With the substitutions 

z A = t 3 U , z 3 = t 3 (1 - ^4) , z 2 = (1 - £ 3 ) *5 , 2i = (1 - £3) (1 - h) (36) 
we finally obtain 



G NP = 2r(3+ r 2 a- ( 2 6 £ ) )r(1 ^ / ^^WM^Il-is)^ [^(tr _3 " 2e 







r(i - 26) 

^(t) = —S\ 2 £3^5 — Ti^t^is — P 3 2 ^^4^5 — P4 ^3^4^4 , U = 1 — t{ 
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In this form the integrand can be fed into the sector decomposition routine for 
"general integrands". The corresponding template file tempi at exbox.m can be 
found in SecDec/general/demos. 



8.3 Phase space parametrisation of a 2 — >• 3 phase space 



The D— dimensional phase space for pi + p 2 — > p 3 + Pa + p$ is given by 



[n 



d u pj 



(2< 



D-l 



5 + {p\ - ml)5 + {pl - ml)5 + {pl - mj)(2n) D 5^ ( Pl + p 2 - £^ 



3=3 



Using 



J d D p 3 5 + {p) - mj) = J d - 1 ^ dE 3 6{E 2 - p 2 - m 2 ) Q{E> 



dn 



D-2 



D-l 

2ix— 

T( D-l\ ' 



dEj dn [ il 2 \p 



\D-3 



\Pj\=y/Ef-^ 



and eliminating p 5 by momentum conservation, one obtains 



/ d$3 = (2ir)™-* \ / dEs ^ dQ °~ 2 ^ D ~ 3 dQ °- 2 b " 4|D ~ 3 

<K(Pi +P2 -P3 -Pif -ml) . (37) 

Having NNLO phase spaces in mind, let us assume that we have two massless 
particles in the final state (which may become unresolved), so m 2 = ml = 0, 
and p^ is the momentum of a massive state, m\ = m 2 , either a single particle 
or a pseudo-state formed by n additional momenta pi in the final state, i.e. 
P5 — Z)r=5P«- I n this case we can parametrise the momenta pi . . .p& as 



Pi 



P2 



^ (i,o^- 2 >,-i) 



p 3 = E 3 (1, 0^ D 4 \ sin 77 sin #i, cos 77 sin #i, cos#i) 

= E 3 (l,0( D - 4 ),n 3 ) 
p 4 = ^4(l,0 (D - 4) ,n 4 ) 



(38) 
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where we choose n 4 such that 



n 4 = 



' 1 

cos 9\ sin 9\ 

V 



\ 



— sin 6*i cos 9\ , 



sin sin # 5 



V 



Thus 



df2gL 2 = d(cos0i) (1 - cos 2 9 1 )^dtt% 



sin (j? sm c/ 2 
cos sin # 2 
cos #2 J 



(</) 



sin sin # 2 



\ 



cos Q\ cos sin # 2 + sin c 4 



d\ cos 9- 

, cos 6 1 ! cos # 2 — sin 6*i cos sin # 2 ^ 



(4) 
D-2 



d(cos9 2 ) (1 — cos 2 9 2 



sm i 



iD-4 



dQ 



D-4 ■ 



(39) 



Due to overall rotational invariance around the beam axis, we can integrate 
out the azimuthal angle r] and use 77 = in eq. fl38|) . leading to n 3 ■ n 4 = cos 6*2. 
We obtain 



^3 = 7 /o ^z^ ^£)-3^D-4^(cosgi) d(cos 2 ) d(cos 0) dE 3 Q(E 3 ) dE i Q(E i ) 



4 (2tt) 2D - 3 



g-4 

sin 2 ^! sin 2 6*2) 2 (sin 2 - 



D-5 
2 



(-E3-E4 

5{m 2 - (p! + p 2 - p 3 - P4) 2 ) 

5(m 2 - ^ +p 2 -Pa -P4) 2 ) = 5(a - m 2 - 2^ (£3 + E 4 ) + 2 £ 3 £ 4 (1 - n 3 ■ rt 4 )) . 
Now we have to choose which variable to eliminate using the 5-constraint. In 
our example, we will eliminate E4, leading to 



(40) 



E d 



s — m 2 — E 3 
2(^i -E 3 {l-n 3 - n A )) 



(41) 



The constraint ®(E±) leads to 
s — m 2 y/s 



77 max 



m 



— j3 where (3 = 1 



(42) 



We substitute 



E 3 = ^Px 3 



E d 



2 H 1 



x 3 



(3x 3 (l - x 2 ) 

cos 9\ = 2 x\ — 1 , cos 9 2 = 2 x 2 — 1 , cos = 2 x 4 — 1 



to obtain 
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d®3 = (2 * _ 3 ^_3^-4 s D - 3 2 D - 8 /3 2D ~ 5 

iD-3 



JJdXi -x{)x 2 (l - x 2 )\ 2 [x 3 (l-x 3 )] J 
i=i 

[x 4 (l-x 4 )]^ [1 -/3x 3 (l -x 2 )] 2 -^ . (43) 
The invariants in this parametrisation are given by 

S 13 = s/3x 3 (l -X X ) 

s 23 = -s(3x 3 x 1 

s/3(l-x 3 ) 



s 3i = (3Kx 3 (1 - x 2 ) , K 



1 - (3x 3 (1 - x 2 ) 
_ 1 -/3(1 -x 2 x 3 ) 
1 - /3x 3 (l -x 2 ) 

Sl4 = -if {t- + x 4 (t + - t~)} = -if Si4 

s 24 = —K |m + — x 4 (u + — m~)| = — K s 24 , 
where 



t± = (yzi(l - x 2 ) ± ^(l - xOj (44) 
= (y/O- ~ Xi)(l - X 2 ) ± y/XiX 2 ^j ■ 
Physical singular limits: 



x 3 — > : 3 soft , x 3 — > 1 : 4 soft , 

xi ->■ 1 : 3 || 1 , X! ->■ : 3 || 2 , x 2 ->■ 1 : 3 || 4 . 

We observe that §u has a line singularity at x 4 = 0, xi = x 2 . We can decouple 
the problem at x 4 = from the line singularity X\ = x 2 by the transforma- 
tion 



r(i-* 4 ) ~ / ^ - \ 

t~ + z 4 (t+ - t~) t- + z 4 (t+ - t-) ' v ; 



leading to 



i i 

j dx 4 [x 4 (1 - x 4 )}^ (Su)" 1 = y dz 4 [z 4 (1 - z^t+t - ]^ t~ + z 4 (t + - t~" 
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1 

/ dZ4 [Z4 (1 — Z4)) 2 — X 2 ) — X 2 (l — Xi) 



n4-D 



(^1(1 - x 2 ) - \fx 2 (l - xi)) 2 + Az 4 \jxi(l - xi)x 2 {l - x 2 ) 
Now we split the ^-integration range at x 2 = x\ and remap to the unit cube: 

1 XI 1 

J dx 2 f(x 1 ,x 2 ) = J dx 2 f(xx,x 2 ) + J dx 2 f(xi,x 2 ) 

xi 



(1) (2) 

where we substitute x 2 = x\ z 2 in (1) and x 2 = xx + (l — Xx) z 2 in (2). Using the 
fact that the contribution from region (2) equals the first one if we transform 
z 2 — > 1 — z 2 and X\ — > 1 — X\ and combining with the original phase space 
given in eq. (Il3l) . we obtain 



rf$ 3 — = fo l D g <KWnj>-4 s D ~ A 2 D - 7 P 2D ^ (46) 

1 

J dxidz 2 dxsdz4[z4 (1 — z±)] xf -3 (1 — X3) D ~ A [1 — /3 £3 (1 — xi + xi^)] 3 D 


xf" 4 ^ 2 D ~ 5 [(1 - Xi)(l - Z 2 )(l - Xi + X X Z 2 )f^ 

4-D 



- Xi)(l - Z 2 ) - VI - X 1 + Xi2: 2 ) 2 + 4 2:4 ^(1 - Xi)(l - Z 2 )(l - Xi + Xi^ 2 ) 



5.^ Troubleshooting 



Below we give possible reasons and solutions for problems which may arise 
during use of the program. 

• The function J 7 is zero: 

verify the on-shell conditions onshell={ . . . } in the file my tempi ate .m where 
you defined the integrand. By default, the external legs have been set to be 
light-like (pf = ssp[i] = 0). If you calculate a massless two-point integral or 
a one-scale 3-point integral, at least one scale must be different from zero 
(e.g. set ssp[l] = — 1 for a two-point function with external momentum pi, 
which amounts to factoring out the overall scale). 

Remember that the program by default replaces pf by ssp[i], {p% + Pj) 2 
by sp[i,j}. (This is done in src/deco/calcFU.m). If symbols different from 
Pi are used for the external momenta, the user has to define their numerical 
values in his template file mytemplate.m in the list onshell. 
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Example: for external vectors called p, q, define numerical values for the 
invariants formed by p and q, e.g. onshell={p 2 — > — l,g 2 — > 0,p*q — > —0.5}. 
Alternatively, you can map to the predefined names for the invariants, e.g. 
onshell={p 2 — > ssp[l],q 2 — > ssp[2],p* q (sp[l,2\ — ssp[l] — ssp[2})/2}. 
This latter solution allows you to leave the invariants symbolic and specify 
numerical values only at the numerical integration stage, by assigning the 
corresponding numerical values in param. input. 

The user can check if the functions J 7 , U and numerator look as expected 
by looking at the file FUN.m in the integralname/ subdirectory. 

• The numerical integration takes very long: 

apart from the fact that this is to be expected for complicated integrands, 
other reasons could be 

• the integrand still contains undefined symbols at the numerical integra- 
tion stage because the numerical values for the constants have not been 
properly defined (e.g. values for which J 7 is not of definite sign, respec- 
tively the general function develops a singularity within the integration 
range). Things to do are: check the function J 7 in the file FUN.m in the 
integralname/ subdirectory (in the loop case); check the log files of the 
numerical integration in the subdiretories 

integralname/polestructure/epstothe [i] , where "polestructure" is of 
the form e.g. 210h0, denoting 2 logarithmic poles and linear, higher 
poles. 

• you chose a very large number of Monte Carlo points and/or a very large 
number of iterations in the input file 

• for functions defined in SecDec/general: verify if there is a singularity 
for Xi — > 1 rather than only for Xi — > and if yes, split this variable at 
1/2 by adding its label to the split list. 

• the results do not appear in an editor window: 

either you did not specify an editor in param . input (last entry) or your sys- 
tem is unable to open the editor window. In this case just look at the result 
file located in the integralname subdirectory (where integralname is the 
name for the calculated integral or graph, specified by you in param. input, 
third item). The result file is called integralname. [point] full . res. 

• you get the message "path for Mathematica not automatically found" : 
Insert the path to Mathematica on your system manually for the variable 
$mathpath in the file perlsrc/mathlaunch.pl. 
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